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ABSTRACT 

Aim of the study. This study examines numerical methods for solving the problems in gas 
dynamics, which are based on an exact or approximate solution to the problem of 
breakdown of an arbitrary discontinuity (the Riemann problem). Results. Comparative 
analysis of finite difference schemes for the Euler equations integration is conducted on 
the basis of an exact or approximate solution to the problem of an arbitrary discontinuity 
breakdown. An approach to the numerical solution of the Euler equations governing 
inviscid compressible gas flow is developed on the basis of the finite volume method and 
finite difference schemes for flow calculation of various degrees of accuracy. Calculation 
results show that monotonic derivative correction provides numerical solution uniformity 
in the breakdown neighborhood. On one hand, it prevents the formation of new 
extremum points, thereby providing monotonicity, but on the other hand, it causes 
smoothing of existing minimums and maximums and accuracy loss. Conclusions. The 
developed numerical calculation method makes it possible to perform high-accuracy 
calculations of flows with strong non-stationary shock and detonation waves and no non- 
Dhvsical solution oscillations on the shock wave front. 
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Introduction 

In some cases, you must be able to calculate numerically supersonic flows 
with strong shock waves (Bulat & Chernyshov, 2016; Bulat & Volkov, 2016, 
Bulat & Upyrev, 2016). Application for solving such problems of traditional 
numerical methods often leads to unphysical solutions (Volkov, 2014; Bulat & 
Volkov, 2015; Bulat et al., 2015). This study examines the accuracy of the 
developed method for the numerical solution of the Euler equations governing 
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inviscid compressible gas flow based on the finite volume method and finite 
difference schemes for flow calculation of various degrees of accuracy. The S.K. 
Godunov (1959) scheme (Einfeldt, 1988), Kolgan scheme (Kolgan, 1972), Roe 
scheme, Harten scheme, and Chakravarthy-Osher scheme (Osher, 1984; Osher 
& Chakravarthy, 1984; Sod, 1978) are used in the calculations; the order of 
accuracy of the schemes varies from 1st to 3rd (Capdeville, 2011). The 
comparison of various finite difference schemes in terms of accuracy and 
efficiency is demonstrated by the calculation of inviscid compressible gas flow in 
a Laval nozzle, with and without a starting shock wave (nozzle shock wave). 

Model problems serve as a testing ground for verifying new methodological 
concepts and assessing the accuracy of the results obtained using software tools 
based on these concepts. An example is checking the developed concepts on the 
projection-evolution method, which is widely used in computational gas 
dynamics (Toro, 2009; Godunov, 1959; Kulikovskii, Pogorelov & Semenov, 2001; 
Roe, 1981). Calculation data provide the means for estimating the monotonicity 
and accuracy of the numerical method and for determining the presence of 
numerical diffusion and non-physical oscillations in the areas of steep gradients 
in the required functions (Donat & Marquina,1996). The problem of breakdown 
of an arbitrary discontinuity (the Riemann problem) has widespread application 
in the finite volume method when it comes to testing computational procedures 
and checking the accuracy of flow calculation schemes (Bulat & Bulat, 2015; 
Bulat et al., 2015). 

To stop the iterative process, the obtained residual level is compared with 
the given degree of accuracy (convergence on machine accuracy is desirable but 
practically unattainable.) Accuracy estimation methods are based either on 
graphical representation of the history of iterative process convergence or on 
theoretical examination of residual behavior. These methods depend on the type 
of convergence (monotonic, oscillating, combined). For this purpose, a hierarchy 
of grids with various decreasing space and time steps is used. 

Grid dependence of the solution can be checked by solving the problem on a 
sequence of grids, where the step decreases by a certain value (by half, for 
example) moving down the grid hierarchy. 

This study solves numerous problems associated with the simulation of 
supersonic flow of an inviscid compressible gas through a Laval nozzle, with and 
without a starting shock wave. The obtained model numerical solutions are then 
compared with the published model solutions, which makes it possible to assess 
the accuracy of the finite difference schemes. Monotonic derivative correction 
provides numerical solution uniformity in the breakdown neighborhood. On one 
hand, it prevents the formation of new extremum points, thereby providing 
monotonicity, but on the other hand, it causes smoothing of existing minimums 
and maximums and accuracy loss. 

The Riemann problem 

The problem of breakdown of an arbitrary discontinuity consists in solving 
the Euler equations in the interval — oo < x < +co under specific initial conditions 
characterized by a constant state in the half-plane -oo < x < 0 (index L) and in 
the half-plane 0 < x< +oo (index R). The initial conditions for the Euler equations 
are as follows (physical variables are used): 
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where is a vector of gas-dynamic variables (density, velocity components, 
total energy per unit of volume). 

A membrane positioned at the point x = 0 separates two gases at different 
pressures, with different densities and velocities. Air with specific heat ratio y = 
1.4 is chosen as the operating medium. 

At the time t = 0, a discontinuity of gas-dynamic parameters is formed at 
the point x = 0, which corresponds to the position of the splitter plate. At the 
time t = 0, the splitter plate is instantly removed. The arbitrary discontinuity 
breaks into several discontinuities, each of which is a shock wave or an 
expansion wave, depending on the initial conditions. Possible solutions contain 
an expansion fan, contact discontinuity, and a shock wave, which split the area 
into four subareas with constant parameter values. The exact solution to the 
problem of breakdown of an arbitrary discontinuity comes down to solving a 
system of nonlinear algebraic equations derived from the conservation law. This 
solution is examined in many research papers (Toro, 2009; Kozhemyakin, 
Omel’chenko & Uskov, 1999). 

Existing solution methods 

In the S.K. Godunov (1959) method gas parameters are approximated by 
piecewise constant distributions on the selected grid such that within each cell, 
these parameters remain constant and equal to the average values over the cell. 
Piecewise constant field evolution over a fairly short time interval, which is 
determined using the exact solution of the Riemann problem in each cell, is used 
to find average cell values on a new layer in time. By repeating this procedure 
step by step, the dynamics of flow variation in time are calculated. 

Piecewise constant and piecewise polynomial distributions of functions in a 
discrete cell, with certain restrictions on the coefficient values of the 
corresponding polynomials, are used to develop highly accurate S.K. Godunov 
(1959) -type numerical methods. One of the difficulties encountered is 
uncertainty in selecting slope values for these distributions. In these cases, the 
breakdown of an arbitrary discontinuity causes it to lose its self-similarity. 
Finding the exact solution to the Riemann problem with random initial 
conditions becomes complicated (for the generalized Riemann problem, in which 
flow parameters on the left and right side of the discontinuity are variables). 
The solution of a non-self-similar problem can be avoided by instead solving a 
problem that is an approximation of the self-similar one. 

The use of piecewise linear approximations requires derivatives to be 
determined in each computational cell. In smooth flows, another approach can 
be applied, where derivatives of each layer in time are approximated by the 
calculated average values. Applying this approach to flows with discontinuities 
results in major errors. In practice, derivatives are not approximated but are 
calculated using data on the nature of the flow. 

The S.K. Godunov (1959) method can be relatively easily generalized to a 
multivariate case. Flow parameters in each grid cell are considered constant (the 
flow in a cell is assumed to be homogeneous). Flow development comes down to 
the interaction of homogeneous flows at cell edges (ignoring flow interaction at 
cell apexes). This approach involves calculating the interaction of two 
homogeneous gas flows, which are initially separated by a certain plane. The 
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solution of a 3D problem comes down to the solution of a one-dimensional 
problem along the normal to the separating plane. 

A previous study (Kulikovskii, Pogorelov & Semenov, 2001) introduces an 
exact solution to this problem. Finding such a solution is crucial because it 
makes it possible to consider the properties of a fairly standard medium when 
an approximation by two-term equations of state is applied. One-parameter 
approximations of the equation of state do not ensure the transmission of 
properties of various media on the left and right side of the discontinuity. For 
some types of equations of state, the solution of the Riemann problem is not the 
only one. The solution has a complex, non-classical structure or disturbs the 
hyperbolicity of a system of equations 

The approaches described in (Roe, 1981; Einfeldt, 1988; Donat & Marquina, 
1996) are among the most widespread methods that are based on approximated 
solutions to the problem of breakdown of an arbitrary discontinuity. In some 
studies (Osher, 1984; Osher & Chakravarthy, 1984), when an arbitrary 
discontinuity breaks down, a shock wave is replaced by a compression wave, 
which results in a system of equations with monotonic solutions. In another 
study (Einfeldt, 1988), breakdown is treated using two Jacobians, one for the left 
node and one for the right node, to construct a single Jacobian from the simple 
waves of the left Jacobian moving to the right and those of the right Jacobian 
moving to the left. A previous study examines only expansion and compression 
waves. The Roe scheme (Roe, 1981) is based on assigning explicit formulas to 
determine linearized values that make up a Jacobian (primitive variables are 
used). However, the Roe scheme has a drawback: it allows the presence of an 
expansion shock wave at the sonic point. This can be avoided by adding extra 
viscosity, thereby modifying eigenvalues in the neighborhood of sonic points. 

The Roe method (Roe, 1981) is based on an exact solution to the Riemann 
problem for a specially linearized system of equations. The solution consists of 
moving discontinuities separated from each other by areas with constant 
variable values. Such a solution is unusual because it preserves the non-linear 
Rankine-Hugoniot relations on a single shock wave and the relations on a single 
tangential discontinuity. The Roe method provides a means for building finite 
difference schemes for conservative hyperbolic sets of equations. 

In the S. Osher (1984) method, an approximate solution to the Riemann 
problem is built for a quasilinear system of equations and is basically a 
combination of Riemann waves only. 

Numerical methods 

Mathematically, the problem of breakdown of an arbitrary discontinuity is a 
Cauchy problem with initial conditions for conservation laws that determine the 
movement of a compressible gas, with the initial distribution of gas parameters 
in the form of piecewise constant functions. 

In the S.K. Godunov (1959) method, piecewise constant distributions of 
functions are used to describe an instantaneous state of a moving medium (flow 
parameters are considered constant within each control volume). Further 
development of flow approximation, which consists of many elementary 
homogeneous flows in time, is determined by the solution of the Riemann 
problem at the edges of control volumes. It is possible to describe the whole 
variety of discontinuity configurations using a relatively simple scheme. 
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Discontinuity velocities and flow parameters in smooth areas between 
discontinuities can be determined by simple mathematics. 

The S.K. Godunov (1959) finite difference scheme is a conservative scheme 
of first-order accuracy. It is widely applied in various problems involving 
numerical simulation of compressible gas dynamics. The S.K. Godunov (1959) 
scheme uses approximate viscosity, which means that artificial viscosity is not 
required for strong discontinuity calculations. During the calculation of weak 
discontinuities such as expansion waves, approximation error becomes 
significant, which manifests itself as strong smearing (the less the Courant 
number, the stronger is the smearing). The restriction on time step follows from 
the condition that waves formed after the discontinuity breakdown at the edge of 
the control volume should not reach its center or, with a less strict restriction, 
its other edge. In the area where strong shock waves interact, the S.K. Godunov 
(1959) method displays numerical damping properties similar to those of the 
application of artificial quadratic viscosity (Kulikovskii, Pogorelov & Semenov, 
2001). The S.K. Godunov (1959) methods based on the solution of the Riemann 
problem by the Roe method have numerical linear viscosity. During the 
simulation of strong shock waves, the Courant number needs to be decreased 
from 1 to 0.2 to ensure stability of calculations. 

The data required for flow calculation in the Godunov method come down to 
using the state that formed after the discontinuity breakdown at the edge of the 
control volume. Solving the Riemann problem is demanding; yet, the Godunov 
method only uses a part of the obtained data. 

Godunov’s numerical methods constructed using the exact solution of the 
Riemann problem provide the means for calculating shock waves of arbitrary 
intensity with the Courant number close to 1 (0.8 < C < 1). 

Numerical methods based on approximate solutions to the Riemann 
problem rely on separate exact elementary solutions to the Riemann problem 
(moving discontinuities in the Roe method or Riemann waves in the S. Osher 
(1984) method). Various conclusions and additions are discussed in previous 
studies. 

Nozzle flow 

Let us consider the flow of inviscid compressible gas in a channel with a 
variable cross-section area (a Laval nozzle). The nozzle section is described by 
the dependence y = [(1 + x2)/n]l/2, where x is in the interval [-0.3.1]. Air (y = 
1.4) is taken as the operating medium. Details on the computational procedure 
have been provided in previous studies (Volkov, 2005; Volkov, 2006). Indices 00 
correspond to the parameters of the flow in the receiver, while the index oo 
corresponds to the parameters of the flow in the surrounding medium 

The nozzle flow mode is determined by the correlation between the pressure 
at the nozzle exit, p2, and the pressure in the tank, pi. In alternative 1 
(pressure drop below critical, p2/pl > 0.528), a flow mode with back pressure is 
implemented, where pressure buildup of the over-expanded flow to the external 
pressure is performed through the nozzle shock wave. In alternative 2 (pressure 
drop is above critical, p2/pl < 0.528), gas is sustainably accelerated from 
subsonic velocity at the point of entry to a certain velocity at the throat, which is 
dependent on the given pressure drop. The gas is then stagnated. 
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At the nozzle entry, the total pressure and total temperature are set. The 
settings for boundary conditions at the nozzle exit depend on the gas flow mode. 
With subsonic flow at the exit, the static pressure is set equal to the external 
pressure. Other parameters are determined through extrapolating variables 
from the internal cells of the computational domain. Supersonic flow requires no 
additional boundary conditions at the exit. 

Calculations are performed on a 100-cell grid. The Courant number is 
assigned the value CFL (Courant-Friedrichs-Lewy) = 0.95. A stationary 
solution to the problem can be achieved using the relaxation method. 

In alternative 1 (p2/pl > 0.528), stagnation pressure (pi = 106 Pa) and 
stagnation temperature (T1 = 300 K) are set at the nozzle entry and static 
pressure (p2 = 8*105 Pa) is set at the nozzle exit. At the entrance to the 
computational domain, the flow is subsonic. At the convergent section of the 
nozzle, the flow is accelerated, reaching sonic velocity at the throat, and it keeps 
moving at supersonic velocity. At the divergent section of the nozzle, a straight 
shock wave is formed, beyond which the flow becomes subsonic. 

In alternative 2, gas is accelerated at the subsonic section and stagnated in 
the supersonic section. 

Calculations are required to compare various flow calculation schemes 
(Yeom & Chang, 2013; Su, 2014). Results of the numerical simulation processed 
in the form of pressure—coordinate x-dependence are shown in the Fig. 1 
(Courant number 0.98) and Fig. 2 (Courant number 0.25). Finite difference 
schemes smear out the discontinuity across 1-2 computational cells, preserving 
the monotonic nature of the solution. 

In the presence of a nozzle shock wave (Fig. 1), the Chakravarthy-Osher 
scheme produces results closest to the exact solution. The Godunov scheme 
requires 0.03479 s of estimated time (12 time steps) to reach a stationary state, 
the Kolgan scheme requires 0.03401 s (12 steps), the Roe scheme requires 
0.01999 s (6 steps), the Harten scheme requires 0.01975 s (5 steps), and the 
Chakravarthy-Osher scheme requires 0.01968 s (5 steps). 

In a supersonic flow mode (Fig. 2), all finite difference schemes produce 
similar results. The Godunov scheme requires 0.007257 s of estimated time (14 
time steps) to reach a stationary state, the Kolgan scheme requires 0.008151 s 
(15 steps), the Roe scheme requires 0.009545 s (18 steps), the Harten scheme 
requires 0.006941 s (12 steps), and the Chakravarthy-Osher scheme requires 
0.006938 s (10 steps). 
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Figure 1. Pressure distribution along the x coordinate with p 00 = 10 6 Pa, T 00 = 300 K, p~ = 
8* 10 5 Pa. 

Icons ■ and □ correspond to the results obtained using the Godunov scheme and the Kolgan 
scheme, respectively (a); ► - using the Roe scheme (b); ◄ - using the Harten scheme (c); 

• - using the Chakravarthy-Osher scheme (d) 
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Figure 2. Pressure distribution along the x coordinate with p 0 o=10 8 Pa, T 00 =300 K, poo=8-10 5 
Pa. 

Conclusions 

It is possible to assess the accuracy and convergence rate of a scheme by 
comparing solutions of model problems with exact solutions. Construction of 
exact model solutions is a crucial element in the general routine of numerical 
algorithm construction. 

Comparative analysis of finite difference schemes for the Euler equations 
integration is conducted on the basis of the exact or approximate solution to the 
problem of an arbitrary discontinuity breakdown. The accuracy and 
effectiveness of various finite difference schemes are demonstrated by 
calculating the inviscid compressible gas flow in a Laval nozzle. The Godunov 
scheme, Kolgan scheme, Roe scheme, Harten scheme, and Chakravarthy-Osher 
scheme are used in the calculations; the order of the schemes varies from 1st to 
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3rd. Calculation results show that monotonic derivative correction provides 
numerical solution uniformity in the breakdown neighborhood. On one hand, it 
prevents the formation of new extremum points, providing monotonicity, but on 
the other hand, it causes smoothing of existing minimums and maximums and 
accuracy loss. 
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